!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! calculate x4, y4 from given nozzle wall function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! wall function is y=a+b*x+c*x^2+d*x^3
! slope of wall is dy/dx=tan(theta)=b+2*c*x+3*d*x^2
! slope of the slope of the wall is d^2 y/dx^2==2*c+6*d*x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   a,b,c,d    wall equation constant
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   all angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine Wall3()
   use VariableDef
   implicit none
   real*8::a,b,c,d
   real*8::da,de,dda
   
!   attachment point
   xa=rtd*sin(thetaa)
   ya=yt+rtd*(1-cos(thetaa))
!   constant
   a=-1*(4*da*xa**3+2*de*xa**3-dda*xa**4-12*da*xe*xa**2+4*dda*xe*xa**3+6*da*xa*xe**2&
      &-3*dda*xa**2*xe**2-6*ya*xa**2+12*xa*xe*ya-6*ya*xe**2)/(6*(xe-xa)**2)
   b=(de*xa**2-2*da*xa*xe+dda*xe*xa**2+da*xe**2-dda*xa*xe**2)/(xe-xa)**2
   c=(2*da*xa-2*de*xa-dda*xa**2+dda*xe**2)/(2*(xe-xa)**2)
   d=(de-da+dda*xa-dda*xe)/(3*(xe-xa)**2)
   
   ye=a+b*xe+c*xe**2+d*xe**3

   if (abs(thetae-thetaa) <= 0.01) then
      c=0.0
      d=0.0
      x4=(a-y2+lp*x2)/(lp-b)
      y4=a+b*x4+c*x4**2+d*x4**3
      theta4=atan(b+2*c*x4+6*d*x4**2)
   else
      x4=(lp-b-sqrt((lp-b)**2-4*c*(a-y2+lp*x2)))/(2*c)
      y4=a+b*x4+c*x4**2
      theta4=atan(b+2*c*x4)
   end if

end subroutine Wall3
